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AN ALGORITHM FOR A GENERAL CLASS OF ROUTING PROBLEMS 


DERIVED FROM HUYGENS' PRINCIPLE 

By Lee M. Avis and George R. Young 
Langley Research Center 

SUMMARY 

If a set of N points or nodes with a nonnegative cost associated with each ordered 
pair is known, it is desired to find a path from one given node to another given node which 
minimizes the cost sum. An algorithm is presented which yields a global minimum solu- 
tion after at most N - I iterations or on a typical large third- gene ration computer, after 
1 hour of computation time for a 10 000- node problem. The rapid- access data storage 
capacity demanded by the algorithm is approximately 3N words for costs read in from 
slow- access storage or 2N words for calculable costs. The time -storage requirements 
of the algorithm compare favorably with those of dynamic programing or any other optimal 
path algorithm known to the authors. When the problem is viewed as a discretized optimal 
control problem, after N - 1 iterations, an optimal control or node transition is estab- 
lished for each of the N nodes or states; thus, the algorithm can be applied to situations 
where there may be errors in the control that necessitate a closed loop control philosophy. 

INTRODUCTION 

The following class of routing problems is considered: If a network of nodes with a 
nonnegative cost associated with each ordered pair of nodes is given, it is desired to find 
a minimum cost path from one given node to another given node. 

A method of finding global optimum paths through a network of nodes is presented as 
a discretized algorithmic form of the Huygens' Principle (ref. 1) which can be stated as 

An optical wave front propagates through space or matter as if each 
point on the front were a source of a new wave front and the envelope 
of these new (secondary) wave fronts in advance of the given (primary) 
front were a subsequent primary wave front. The line joining a sec- 
ondary source and the envelope of secondary wave fronts in the shortest 
optical distance approximates a primary optical ray; in general, the 
nearer the envelope to the secondary source, the better the approximation. 




The algorithm presented herein constructs (exact) optimal paths through the network of 
nodes as the Huygens T Principle constructs light rays through space or matter, where 
optical distances and propagating wave fronts are replaced by, respectively, cost dis- 
tances and propagating cost contours. 

The algorithm is presented and developed as a means of solving the fixed end points 
optimal path problem stated in the section "Statement of the Problem"; extensions to vari- 
able ends points and optimal control problems are discussed in the section "Applications 
to Optimal Control Problems With Multinode Goals and Origins." 

The Huygens' Principle foundation of the current algorithm should lead to much fur- 
ther generalizations and extensions of the technique developed here than that of a similar 
algorithm developed by E. W. Dijkstra (ref. 2). 

SYMBOLS 

a,b computer times in microseconds for certain operations 

c(I 0 , . . .,j) cost of path (Iq, . . .,j) 

(Iq, . . .J) path from node Iq to node j 

(Iq — j) optimal path from node Iq to node j 

(i,j) direct path from node i to node j 

I 0 ,I f starting node and goal node, respectively, defining the end points of the optimal 

path 

I_ starting node defined on page 9 

(k + l)th node to be occupied; defined on page 10 
j k unoccupied node for which c(I 0 , . . .,j) is a minimum at kth step 

cost matrix element; the cost of the direct path from node i to node j 
N number of nodes in problem 

indices of refraction 
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q any node on path of equation (2) 


r,s,t,j running indices 


( s k) 


(k + l)th set of occupied nodes, {Jo^l’^2’ * * *>*^3 


T k cost of an optimal path from node Ig (or Jg) to node 


AT 




J 


T k+ i - T^ at the (k + l)th iteration of algorithm, which is obtained opera- 
tionally as minimum of r values 

next to last node of optimal path from node Ig to node i for any occupied 
node i; for an unoccupied node i, is a current estimate of the next to 
last node of optimal path from node Ig to node i 

at (k + l)th iteration, the cost distance remaining to reach (occupy) node j 
along the direct path (J r ,j) where J r is an occupied node (defined after 
eq. (4)) 

at (k + l)th iteration, the minimum over all occupied nodes J r of 

(kl 

for node j unoccupied; for node j occupied, t > ' = <*> 
variable of algorithm in step (6) on page 7 


(k) 

J r ,J 


STATEMENT OF THE PROBLEM 


A set of N nodes, numbered arbitrarily from 1 to N, together with a cost matrix 
K- j is known where i,j = 1, . . .,N and = 0; find a sequence of nodes 


f-1 


,1^) which minimizes 


i=0 


K 


Vi 


i^l 


when Iq and I f are known but not the number 


of nodes in the path f + 1. That is, a minimum cost path from Ig to If is sought. A 
sequence of nodes is called the ’’path" through those nodes and the cost matrix defines the 
cost of all paths. 


EXAMPLE PROBLEM 


A simple four -node example problem is solved to illustrate the principles of the 
algorithm and its kinship to Huygens’ Principle. Consider the four nodes at the vertices 
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of the parallelogram of sketch (a). The costs associated with the sides and diagonals of 
the parallelogram are indicated in sketch (a). Suppose an optimal path from node 1 to 
node 4 is desired. The cost per unit path length is n 2 to the right of (but excluding) the 
diagonal (3,2) and n x elsewhere. If n i and n 2 are interpreted as indices of refrac- 
tion, then the light ray path from node 1 to node 4 of minimum optical path length (or 
equivalently, minimum travel time) is sought when the ray is confined to the sides and 
diagonals of the parallelogram. A discretized simulation of the refraction of light rays 
passing from one medium to another of different refractive index has been defined. 


n l I n 2 



Let n t = 1 and n 2 = 3. The matrix of costs Kjj associated with the ordered node 
pairs or links is then 


0 



25.30 


6 7.21 25.30 

0 4 21.63 

4 0 18 

21.63 18 0 


Consider a wave front originating at node 1 and propagating at a speed inversely propor- 
tional to the cost per unit link length; that is, the wave front is a contour of constant cost. 
In order to indicate that node 1 has been encountered by the wave front or "occupied," 
is set equal to a very large number (written °°), The shortest cost distances to the 
unoccupied nodes 2,3, and 4 from the wave front along direct links to these nodes are 
specified by, respectively, = 6, = 7.21, and = 25.30. To complete the descrip- 

tion of the wave front, record the direct link associated with each t. by uj = i where j 
is the node and the direct link is (i , j) . Thus, at the beginning, u i - u 2 = = u 4 = 1. 
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The next node reached (occupied) by the wave front is node 2, since r 2 is the mini- 
mum t. in reaching node 2 , the wave front has advanced or six cost units, and so 
the r values of the unoccupied nodes are decremented by 6 . Then is set equal to 
00 to flag node 2 as occupied. Thus the r values become 

Tj = °o r 2 = °° r 3 " r 4 = 19.30 

The newly occupied node 2 now becomes the source of a secondary wave front that is 
consistent with Huygens' Principle. As^ shown in sketch (b) and by the K matrix, the : 
wave front originating at node 2 (now a single point at node 2 ) is more cost distant from 
the unoccupied nodes 3 and 4 than is the wave front originating at node 1; therefore, 
is not changed to nor is u^ changed to 2 , and r 3 is not changed to nor is 

U3 changed to 2. 


u 3 = l u 4 = l 

r r 3 = 1.21 f' Tj = 19.30 



T 1 = °° T 2 = °° 

Uj = 1 u 2 = 1 

Sketch (b).- Node 2 is occupied. 


The possibility that the optimal path from node 1 to any unoccupied node passes 
through node 2 has thus been eliminated, and the secondary wave originating at node 2 can 
be neglected. The Huygens' Principle has a provision similar to that exercised here: 
that all portions of the secondary waves which are behind the advancing wave envelope are 
neglected. 

At this stage, is the minimum t; thus, node 3 is the next node reached by the 
(composite) wave front. By proceeding as before, the new is set equal to r 4 ■ ^3 = 
19.30 - 1.21 = 18.09 and then t 3 is set equal to Now node 3 becomes a new wave 
source. Since this third wave is slightly "closer" to node 4 than is the old composite 
wave, is changed from 18.09 to = 18 and u^ is changed from 1 to 3. See 
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sketch (c). Finally, t 4 being the minimum r, node 4 is the next node to be reached by 
the composite wave. An optimal path from node 1 to node 4 is defined by 

U 4 = 3 


or path (1,3,4). See sketch (d). The optimal paths from node 1 to node 2 and to node 3 
have also been found; u 2 = 1 means that path (1,2) is optimal (once r 2 becomes «>) 
and ug = 1 means that path (1,3) is optimal. The u values define an optimal transition 
to each occupied node for a given starting node. If the problem is worked in reverse from 
node 4 to node 1 (reversing the indices of K if K is asymmetric), the u values would 
define an optimal transition from each occupied node for a given goal node. 


u 3 = 1 u 4 ■ 3 



Sketch (c).- Node 3 is occupied. 


= 3 


Wave 
" front 



Sketch (d).- Node 4 is reached. 


6 



The Huygens’ analogue to the u values are those vectors (rays) connecting each 
secondary source to the envelope of secondary waves in the shortest optical distance. 
The algorithm replaces optical distance with cost distance and replaces the rays with 
pointers to the sources, that is, the u values. The envelope of secondary waves is 
analogous to the composite wave or the cost contour of the algorithm. A distinction 
between Huygens’ Principle and the algorithm is that light, according to Huygens’ 
Principle, propagates so that obstructions cast shadows, whereas the cost contour con- 
structed by the algorithm spreads throughout all available "space.” (The analogue to 
an obstruction is a set of nodes which can be reached from nodes outside the set only 
through links of very high cost.) 


DESCRIPTION OF ALGORITHM 


The algorithm is now presented in a form convenient for calculation. The elements 
of K, the cost matrix, as defined previously, can either be calculated as needed or pre- 
computed and stored in auxiliary storage and read in as needed. Let Iq be the starting 
node and If be the goal node and N be the number of nodes. The algorithm follows. 


Step 

Operation 

Comment 

(1) 

W V° o 

Initialization of iteration loop. 


(2) 

% " 

Flags 1^ as occupied. 

. 

(3) 

AT (= 7 j 0 ) = min jr.J (j * b ■ ■ -,N) 

Begin iteration loop. 


(4) 

•If jo - If. H° to step (10) 

Stop iteration when goal node If 

is occupied. 

(5) 

T. = 

Flags Jq -as occupied. 


(6) 

^; = K i a = i, . • • >n) 
1 2 o’’ 

t 

If K elements are not calculated, read j^th row of K into t\ 1 
Skip this operation if K elements are calculated. 

(V) 

T. - T - AT (All T. t on) 

J J k J / 

• 


(8) 

t. - min < t. , t! > ( All t. t <*) 

J 1 i )J \ ] J 

II K elements are calculated, the rj array is not needed and 
can be replaced by K- . to reduce data storage by N words, 

(9) 

u- =- Jq J for all j having t. = rjj. Go to step (3), 



(10) 

Output j Q ; j = j Q ; Cost - 0.0 

Initialization ol output loop. 


(11) 


Output loop (steps (11) to (15)). Yields an optimal path from I Q 



to L (in reverse order), the total path, cost, and the incremental 

(12) 

Cost - Cost + 

path cost. 


(13) 

Output i, cost 



(H) 

11' i = Iq, exit 



(15) 

j = 1. Go to step (11). 

__s 

/ 
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APPLICATIONS TO OPTIMAL CONTROL PROBLEMS WITH 
MULTINODE GOALS AND ORIGINS 


The optimal path problem is cast in the form of an optimal control problem by 
identifying the nodes as states of a system which is to be brought from an initial state to 
a specified final state through a sequence of states minimizing the sum of the state transi- 
tion costs. Each state transition cost is assumed to be a function only of that state tran- 
sition; this assumption amounts to a definition of the concept of "states, 1 * in that the cost 
matrix defines the system. The states or nodes might represent a discretized approxi- 
mation to a continuum state space. Disallowed state transitions are assessed a high cost 
level to prevent, or at least flag in the output, the "forbidden" transitions. The "forbidden" 
flag is set less than the "occupied node" flag (°°) to avoid premature assignment of nodes 
to the occupied status. 

As pointed out previously, if the optimal path problem is worked in reverse (with 
the goal node as the initially occupied node) with a transposed K matrix, then the algo- 
rithm establishes an optimal transition Uj from each node 3 which becomes occupied. 
The u values thus define an optimal control, or an optimal way to proceed, from each 
occupied node, when a specified goal node is given. Since one additional node is occupied 
upon each iteration, N - 1 iterations are required to generate the complete field of 
optimal controls over the state space. The complete field of optimal controls defines the 
optimal procedure even after a control error, which allows closed loop control with state 
information feedback. 

After N - 1 iterations, all nodes are occupied, including those "isolated" nodes, 
if any, from which the goal node can be reached only by a series of transitions including 
at least one forbidden transition. In this situation, the u values contain some forbidden 
transitions, which, however, are flagged in the output by the associated high cost. In some 
cases, it is convenient and computationally advantageous to remove isolated nodes from 
the state space beforehand. 

The modifications to the algorithm for generating the complete optimal control field 
are 

(1) Replace the K matrix with its transpose [k t ] 

(2) Wherever Iq occurs, replace it with I f , the goal node 

(3) Initialize an iteration number k = 0 before the third step of the algorithm 

(4) Replace the test at step (4) by an operation incrementing k by 1 and then an 

operation testing k for the value N - 1, a successful test triggering a shift 
to the output 
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(5) Replace the output loop (steps (10) to (15)) by a routine which outputs the u- 
array and the associated costs . 

Once the complete optimal control field is generated, the optimal path(s) from any 
node(s) to the goal node If can be output as the sequence (s) of optimal controls beginning 
at the desired starting node(s) I s ; the sequence(s) are of the form ^I s ,Uj , . . . ,1-fj . 

If only the optimal paths from each of a set of starting nodes to a goal node and the 
associated costs are desired (for example, all optimal paths originating within a neighbor- 
hood of a nominal starting node), then N - 1 iterations of the algorithm, in general, are 
not necessary. The modifications to the algorithm to allow less than N - 1 iterations 
are as follows: 


Steps (1) to (3) are the same as those for complete optimal control field 

Step (4) Replace the test at step (4) by an operation incrementing k by 1 when the 
newly occupied node jg is one of the desired starting nodes and perform 
an operation testing k for the number (quantity) of desired starting 
nodes, and then switch to output when k equals the number of starting 
nodes 


Step (5) Replace the output loop, steps (10) to (15), by a routine which outputs the 

sequences of optimal controls ^I g ,Uj s ,u u ^ , . . .,Ifj for each desired start- 


ing node I s and the associated costs 

A generalization of the problem of path optimization with fixed end points, optimiza- 
tion over all paths with starting and goal nodes each in specified respective subsets of the 
set of N nodes, is solvable in one application of the algorithm. A simple procedure is 
to set all K^- =0 for both i and j in the goal subset, transpose the K matrix, set 
Iq equal to any node in the goal subset, change the stopping rule test (step (4)) from 
"j 0 “ If ?" to ,r j 0 e {starting subset}?”, and change the exit rule test in the output loop 
(step (14)) from ”i = Iq?” to M i € {goal subset}?”. 


THEORETICAL BASIS OF ALGORITHM 


In this section, theorems are proven which provide the basis for the algorithm. 

As demonstrated by the example problem (and shown in the following) , the algorithm 
places each node in the occupied status, in sequence, in increasing order of the cost of 
optimal paths from the starting node I 0 to each of the nodes. It is first reasoned that 
closed loops need not be considered in seeking an optimal path; thus, only a finite number 
of paths need to be considered. From this finite number of path costs, a unique minimum 
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can be selected. It is proven that there is an optimal path from Iq to an unoccupied 
node j such that (a) all nodes in the path except j are occupied and (b) the cost of the 
path is not less than the cost of any optimal path from Iq to any occupied node and not 
greater than the cost of any path from Iq to any unoccupied node. The procedure for 
computing the cost distance to node j is then demonstrated. This procedure, applied 
recursively, yields a global optimum path. 

Assume that a ranking (with ties permitted) exists for the costs of all paths orig- 
inating at Iq. The existence of a unique minimum path cost over all paths from Iq to 
any given node is shown to result from the nonnegativity of the elements of the cost matrix 
K by the following statement. 

For any path from Iq to a given node j containing any node q more than once 
(that is having a closed loop), (I 0 ,n 1 ,n 2 , . . .,q,m 1 ,m 2 , . . .,q, . . .,j), there is a path 
from Iq to j of no greater cost obtained by deleting from the node sequence all nodes 
following the first occurrence of q up to and including the last occurrence of q; for the 
nonnegativity of the cost elements implies that the cost of the sequence of nodes deleted is 
nonnegative. The number of paths from Iq to j with no multiply occurring nodes is 
finite and, by assumption, their costs can be ranked; therefore, a path from Iq to j of 
unique minimum cost over all paths from Iq to j exists, this being the minimum over 
the paths from Iq to j with no multiply occurring nodes. 

Let the minima over the costs of all paths from I Q to each of the N nodes be 
sequenced in a monotonically nondecreasing' order: T^T-^Tg, . . ,,T N _ 1 . For each T k , 

associate a node J k so that the cost of an optimal path from Iq to J k is T k for 
k = 0,1,2, . . .,N - 1. (The sequence of J k values is the sequence in which each node is 
placed in the occupied status. This sequence is not necessarily unique because of the 

possibility of ties, which are arbitrarily ordered.) Let^S^be the set of all J r for 
r = 0,1,2, . . .,k; {S k } is to be identified as the (k + l)th set of occupied nodes. For sim- 
plicity, the condition that Jq = I 0 is chosen without loss of generality. 

The following theorem demonstrates the constructive nature of the algorithm. 

Theorem : There is an optimal path from Iq to some j £ { s k} 

of cost, T k+ ^, such that all nodes in the path except j belong to 

(Sk>, f °r k = 0,1,2, . . ,,N - 2 . 

Proof: By definition of (S k }, T k , and J r for r = 0,1, . . ,,k, 

T k+1 = min (c(I 0 -j)} (1) 


10 



where (Iq - j) is an optimal path from Iq to j and c(Io — j) is the associated cost. 
Let be a j for which c(Io — j) is minimized. Then, 

T k + 1 - <*0 - 3k) (2) 

Let q be any node on the path of equation (2). Then, 

T k+1 = c(I 0> • • -A> • • -’3k) 

If (Iq, . . . ,q) or (q, . . .,j^) is not optimal, then (Iq, . , . ,q, . . .,j^) is not optimal. By 
contradiction, 

T k+1 = c Oo - I - ik) 

By the nonnegativity of Kjj , 


c(I 0 - q) s e(I 0 - q - J k ) = T k+1 

Suppose q i (S k ). Then q is a candidate j in equation (1) and 


(3) 


c(I 0 " <1> S T k+1 


Therefore by equation (3) 


c(I 0 - q) = T k+1 


(Note the incidental result that c(q — j^) = 0.)- 

The algorithm determines the path -defining u values indirectly by means of recur- 
sive minimization over the incremental costs (the r values) of completing optimal paths 
to unoccupied nodes. The essential features of the optimal path determination through 
recursive minimization are set forth in the following theorem; the corollary to the theorem 
allows simplification of the computation and reduction of data storage. - T^ is the 

AT of the algorithm at the (k + l)th iteration.) 


Theorem: 


T 


k+1 


- T k = 


min 

re{0,l, . . .,k) 

^{ s k} 




(4) 
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where 


a ( k ) _ a (k~l) - ( t. - Tu 
a j r ,j J r ,j \ k k -v 


(r e{0,l (k- 1)}; j i {S k }) 




( r = k; j^{S k >) 


Furthermore, there is an optimal path from Iq to J k+1 of the 
form fl 0 — J R , J k+1 ) where J R and J k+ ^ are the J r and j, 

' ' (k\ 

respectively, which minimize a > ' .. 

u ^ 


Proof: Equation (1) in the proof of the first theorem is 


T k+1 = min ( c(I 0 - j)} 

j^{Sk> 

which, by the first theorem and the definition of the J r , can be written as 

T k+1 = min { c do ~ J rd)} 
re{0,l-, . . .,k} 

i/{Sk} 

By definition, 

c(I 0 J"r) = T r 

Thus, equation (5) can be written as 


(5) 


T k+1 - min ^ c ^r>j^} 

re{0,l, . , ,,k) 

j/{Sk} 
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By subtracting TV and substituting K, . for c(J r ,j), 
K J r 


k+1 


- T k = 


From the definition of 


min 

re{0,l, . . .,k} 

jX{S k > 

•, it follows that 
d r> J 




( 6 ) 


4 k) . = oW . 

J r d J r ?j 


- (T k * T r ) = K 


J r’i 



Therefore, equation (6) can be written as 


€ <0,1, . . .,k}; j i {S k }) 


T k + 1 • T k = r mil 
re{o,l, . 

i/{s k ) 

which establishes the first part of the theorem (eq. (4).) Let J R and J' be the J r and j, 
respectively, which minimize . in equation (4). Then, by equations (5) and (4), 



T k+1 = <*0 “* j R^ jT ) 


By equation (1), the path (Io — Jr,J ) is optimal; thus, by definition of the J r , one can set 
J k+1 = J * Therefore, (Io - JR,Jk+l) is optimal. 


Corollary; 


T k+1 ' T k = min 




(f 
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where is defined by 




Proof: The proof is inductive on k. 


Suppose 



min 

te{o,l, . . ,,k} 



j 4 (S k ) and some k e {0,1, . . ,,(N - 3)}j. Then, by equation (4) and the hypothesis, 


T k + l' T k= min 
By equation (4), 




- K J r ,j. 


( r = k + 1; j 4 { s k+i}) 


Thus, 


re 


. min . ft j 5 = min (t fe - 1 

(0,1, . . ,,(k+l)} ^ r J L 1 


T k + 1 - T k))> K J k+1 ,j 


where t e (p,l, . . . ,k}, j or by the hypothesis 


n~ 


min 


, ?} = min B k) 1 - ( T k + i ■ T k 

re{0,l, . . .,(k+l)} ^ r J 


K I , 

L k+1’L| 


- T (k+1) 
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Since 


r 


r (k) = 


min 


° { 3 


T k + 1- T k= 


j/{ s k} 


t€{0,l, . . .,k> ^ v V ) 


r (k+l) _ 
j 


min 


re{0,l } . . .,(k+l)} 


ft?} 




for all j £ and for any k f (0,1, . . .,(N - 3)}, in order to complete the inductive 

proof, it will be shown that for k = 0, 


r9$ = min <fcri k ^ 

te{0, . . „k} V t>J 


All j i {S k } 


For k = 0 by definition, 


T (k) _ K 

i _K V 


(ah j i {s k } 


min i a j k) jl = t7 j° ) j = K j i 

u{o, . . „k} J o J J o- J 


Therefore, 


T (k) = 


nun < a- 


(k). 


J te{0, . . .,k} ^ 1 


Note that the algorithm defines the optimal paths by updating the uj when the Tj 
are updated by 




1 

rj k ^ = min<j 

f' 1 * - ( T k - T k-i)j- 

c* 
*“ J ‘ 1 



■ -J 


15 



whereupon if Kj . < t!^ 1 ^ - ^T k - T k _ ^ , then u- is set to J k . Since the proof of the 


corollary shows that 


T ( k ) = 


t£{0, . . .,k) k 1 


(ke{0, . . . ,(N - 2)}; j 4 {S k }j 


then, when the minimum r( k ^ is found, m is the J. which minimizes al k k Thus, 

J J L d t’l 

by the second theorem, there is an optimal path from Iq to J^+i of the form 

(Iq - Uj ,J k+1 j where J k+1 is the j which minimizes rj k \ 


COMPUTATIONAL CONSIDERATIONS 

The rapid- access data storage demanded by the algorithm is approximately 3N words 
(for the uj, Tj, and rjj for cost elements read in from slow-access storage. (The rapid- 
access data storage can be reduced by up to N words by reading at step (6) partial rows 
of the K matrix into an abbreviated r 1 array, working with this array down through 
step (9) with appropriate changes to running index limits, and cycling back through step (6) 
for the next r' array until the j Q th row of K is exhausted. The additional computation 
time that would result might be acceptable if core storage were a serious constraint.) 

For cost elements calculated as needed, the rapid-access data storage is approximately 
2N (for the uj and the Tj). Now consider computation time, or more precisely, resi- 
dence time in core memory. 

The numbers of operations performed for the worst case of N - 1 iterations are 


as follows: 



Operation 

Total times performed 

Breakdown 

Comparisons 

2N 2 

Step (7), N 2 
Step (8), 1/2 N 2 
Step (3), 1/2 N 2 

Additions or subtractions 

1/2 N 2 

Step (7), 1/2 N 2 

Substitutions 

2N 2 

Step (3), N 2 
Step (8), 1/2 N 2 
Step (9), 1/2 N 2 
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Operation 

Total times performed 

Breakdown 

Calculations of r or 

1/2 N 2 

Step (8), 1/2 N 2 

Row read-in of 7,7* from auxiliary storage 

2N 

Step (6), 2N 


For those particular operations where an exact count of frequency could not be made, the 
worst case frequency is taken. 

For large third-generation machines, the computation times of each of the listed 
operations except the last two are in the approximate range 1 to 10 microseconds. The 
time used in a row read- in of r ,r' or a read- in of a single element is of the order of 
100 milliseconds for disk storage systems. Thus, for costs read- in rather than calculated, 
the computation time is approximately 

Time % 4.5 x 10“ ® aN^ + 0.2N seconds 

where a is a characteristic computation time in microseconds for these listed opera- 
tions (excluding the last two), and 1 ~a ~ 10. For the case of computed costs, the com- 
putation time is approximately 

Time ^(4.5a + b/2) x 10" 6 N 2 seconds 

where b is the computation time in microseconds for one 7 computation. 

COMPARISON WITH OTHER ALGORITHMS 

The data core storage and maximum numbers of operations for the proposed algor- 
ithm, a dynamic programing algorithm (ref. 3), and an algorithm due to Ford and Fulkerson 
(ref. 4) are presented in the following table. The term ’’elementary operations” refers to 
replacement, addition, subtraction, multiplication, division, or simple one-branch or two- 
branch logical operations. 

The Ford and Fulkerson algorithm is limited in practice to problems with no more 
than a few hundred nodes because of the rapid increase in the number of operations with 
increasing numbers of nodes. However, this algorithm works under a less restrictive 
assumption than that of nonnegative costs; the cost of any closed loop is nonnegative. It 
can be seen from the table that the proposed algorithm is comparable in storage to the 
Ford and Fulkerson algorithm but requires less computation time by a factor of the order 
of N 2 . It can also be seen from the table that the general routing problem with nonnega- 
tive costs is solvable by the proposed algorithm in, for worst cases, less computer time 
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by nearly a factor of N and with somewhat less core storage than the dynamic program- 
ing algorithm. This advantage can be exploited to solve general routing problems with 
an order of magnitude more nodes than those solvable by dynamic programing in the same 

run time. 


— - 

Algorithm 

Source of cost 
elements 

Data 

core 

storage 

Number of 
elementary 
operations 

Number of 
cost element 
computations 

Number of 
accesses of 
auxiliary storage 

Dynamic 

programing 

Computed 





Dynamic 

programing 

Read in from 

auxiliary storage 

=4N 

*3N 3 


~N 2 

Dynamic 

programing 

Stored in core 

~N 2 

_ 

^3N 3 



Ford and 
Fulkerson 

Computed 

*2N 

—1/8 N 4 

*1/8 N 4 




Ford and 
Fulkerson 

Read in from 

auxiliary storage 

~3N 

—1/8 N 4 

— 

—1/8 N 3 

— — 

Ford and 
Fulkerson 

Stored in core 

~N 2 

— 1 / 8 N 4 



Proposed 

Computed 

*2N 

*9/2 N 2 

—1/2 N 2 


Proposed 

Read in from 

auxiliary storage 

==3N 

=9/2 N 2 



*2N 


CONCLUDING REMARKS 

Huygens’ Principle has been used as the motivation to develop an algorithm for 
solving a general class of routing problems. If a nonnegative cost for each ordered paii 
of a set of N nodes is known, the algorithm can find the path of minimum cost sum from 
one given node to another given node in at most N - 1 iterations. By viewing the piob- 
lem as a discretized optimal control problem, after N - 1 iterations, an optimal conti ol 
or node transition has been established for each of the N nodes or states; thus, the 
algorithm can be applied to situations where there may be errors in the control that would 
necessitate a closed loop control philosophy. The algorithm is compared in both data 
core storage and maximum number of operations with a dynamic programing algorithm 
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and with an algorithm of Ford and Fulkerson. The algorithm is found to compare favor- 
ably in both core storage and maximum operations for a large number of nodes. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., February 25, 1974. 
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